PR Review: Implement the Getis-Ord Gi* statistic in hotspots()

Blockers (must fix before merge)

None.

Suggestions (should fix, not blocking)

• focal.py:1296 _gistar_zscore does its arithmetic in float32 (whatever the convolution returned) and casts at the end. For large n and large weights, n * sq_weight_sum - weight_sum^2 is a difference of large numbers and can lose precision in float32. The hand-reference tests run in float64 and pass at rtol 1e-3 to 1e-4, so the tolerance hides it. Accumulating the variance term in float64 before the sqrt would be safer. Not blocking: the z thresholds (1.65/1.96/2.58) are coarse, so float32 rounding rarely flips a band.

Nits (optional improvements)

• focal.py:1305 _gistar_convolutions_numpy is numpy-only and the tests import it directly, while the dask path inlines the same three convolutions in _hotspots_chunk. If those two spellings drift, the tests only catch the numpy one.
• The dask+numpy and dask+cupy paths each compute n/mean/std with three reductions in one da.compute(...). It is a single pass, but the mask .sum() and nanmean/nanstd rebuild the NaN mask three times. Minor; not worth changing.

What looks good

• Gi* z-scores checked against an independent hand-computed reference on a known window, for binary and weighted kernels, not just the classification banding.
• NaN cells drop out of the neighborhood sums and the global stats via a validity mask convolved with the kernel and the squared kernel. That is the correct treatment.
• All four backends share the same _gistar_zscore math through the existing dispatch. numpy, cupy, dask+numpy, and dask+cupy all pass on a real GPU here.
• n < 2 and zero-std guards have clear errors.
• The numpy/dask edge divergence is the documented boundary='nan' artifact; the tests compare the interior, which matches the existing dask+cupy hotspots test.

Checklist

• Algorithm matches the reference formula
• All backends agree on the interior
• NaN handling is correct
• Edge cases covered (NaN, constant surface, n<2, tiny raster)
• Dask chunk boundaries handled (depth = kernel radius)
• No premature materialization beyond the documented scalar pass
• Benchmark: existing hotspots benchmark still applies; no new function
• README feature matrix: no change needed
• Docstrings present and accurate (formula added)

Follow-up review (after b36e43f)

Disposition of the first-pass findings:

• Suggestion (float32 variance precision, focal.py:1296): fixed in b36e43f. The weight sums and numerator are promoted to float64 before the sqrt, then the z-score is cast back to float32. n * W2 - W^2 no longer loses precision for large rasters or weights. Banding is unchanged on the existing tests; interior z-scores now match the float64 reference to ~2e-7 (the residual is from the float32 convolution inputs, not the variance term).
• Nit (numpy-only _gistar_convolutions_numpy vs the inlined dask convolutions): dismissed. test_hotspots_gistar_dask_matches_numpy_full runs the dask inline path against the numpy reference on the full interior, so drift between the two spellings would fail a test.
• Nit (three reductions rebuild the NaN mask): dismissed, as flagged in the first pass. It is a single da.compute pass and the cost is negligible.

Re-ran test_focal, test_emerging_hotspots, test_validation, test_dask_laziness: 296 passed, including the GPU hotspots paths on a real device. No new findings.